REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No  0704^JSS 


3.  !^~or  .  7  ;^:.  '^'•^■"5  'O'  '«lof o,  t^.,  o^ar-  t:  'V^.no.jo  tO'--J-t^  >r^4-a.'y,  t-.,  O^'OO*  r...".*!,  o-  .t",.  o' 

'■  •  ^fOva.Q*^  »-:x^T  to  ?C<-J:M>  W4vA«>v;to^  CK  20V31 


’  UNIT  fi„v,  to/^no  I  2.  report  0*Tt  1  3.  REPORT  TYPE  *N0  OAT£S  C0VER£0 - 

_ _ _ _ _ June  26,  1995  Final  Report,  5/1/91  -  4/30/9^ 

4.  TITLE  ANOSuiriTLf  - - - - - - - - - — - - 

S.  FUNDING  NUMBERS 

Multi-D linens ional  High  Order  Non— Oscillatory  Numerical 

Methods  for  Discontinuous  Problems  in  Parallel  Structure  DAAL03-91-G-0123 
6.  AUTHOR(S)  - - - - - - 


Chi-Wang  Shu 


7.  PERFORMING  ORGANIZATION  NAM£(S)  AND  ADDRESSEES) 
Division  of  Applied  Mathematics 
Brown  University 
Providence,  RI  02912 


9.  SPONSORING /MONITORING  AGENCY  NAM£(S)  AND  ADDRESSEES) 
U.S.  Army  Research  Office 
P.  0.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


jwms  I 


10.  SPONSORING  -  MONITORING 
AGENCY  REPORT  NUMBER 


Tl-  SUPPLEMENTARY  NOTES  ’  - - - - * - — - 

The  views,  opinions  and/or  findings  contained  In  this  report  are  those  of  the 

author (s)  and  should  not  be  construed  as  an  official  Departaent  of  the  Army 

position,  policy,  or  decision,  unless  so  designated  by  other  documentation. 

OISTRiaUTlOW  y  AVA(LA>IUTY  STAT£MEWT  '  1  ^ 

»i*iewtNi  I  12h.  04STRItUTION  COOC 


Approved  for  public  release;  distribution  unlimited. 


13.  ABSTRACT  ■>rin  .  ■  i 

In  this  project  we  have  studied  high  order  finite  difference,  finite  element, 
and  spectral  methods  for  the  numerical  solution  of  discontinuous  or  high  gradient 
solutions.  Theoretical  study  of  algorithm  development,  stability,  accuracy  and 
convergence  analysis,  as  well  as  applications  of  the  algorithms  to  computational 
fluid  dynamics  (both  compressible  and  incompressible)  and  semiconductor  device 
simulation  are  performed.  Parallel  implementation  Issues  of  these  high  order 
methods  are  also  studied. 


19951026  036 


Wig  QUALITY  INSPECTED  8 


>«.  SUtJECT  TERMS 


IS.  NUMStR  OF  FACES 

8 

1*.  FWCF  COO£ 


'7.  SCCURITY  a_ASSlFK"ATi<-i,.  (  '• — ~  _ 

1  20.  LIMITATION  of  abstract 


OF  RfFORT 

unclassified 

NSN  75aO.Ol. 260-5500 ' 


OF  THIS  PAGt 

UNCLA.SSIFIED 


OF  abstract 

UNCLASSIFIED 


Stand'ifd  Form  29S  ERev  2*89) 

»'r>cr.o<c  t>v  anSj  Sta 


Final  Report  of  ARO  Grant  DAAL03-91-G-0r23 
Multi-Dimensional  High  Order  Non-Oscillatory  Numerical 
Methods  for  Discontinuous  Problems  in  Parallel  Structure 


Chi-Wang  Shu 

Division  of  Applied  Mathematics 
Brown  University 
Providence,  RI  02912 
E-mail;  shu@cfm.brown.edu 

May  1,  1991  to  April  30,  1995 


1.  Abstract 

In  this  project  we  have  studied  high  order  finite  difference,  finite  element,  and  spec¬ 
tral  methods  for  the  numerical  solution  of  discontinuous  or  high  gradient  solutions. 
Theoretical  study  of  algorithm  development,  stability,  accuracy  and  convergence  anal¬ 
ysis,  as  well  as  applications  of  the  algorithms  to  computational  fluid  dynamics  (both 
compressible  and  incompressible)  and  semiconductor  device  simulation  are  performed. 
Parallel  implementation  issues  of  these  high  order  methods  are  also  studied. 


2.  Statement  of  the  Problem  Studied 


The  problems  studied  in  this  project  involve  numerical  solutions  of  partial  dif¬ 
ferential  equations,  mainly  hyperbolic  type  or  convection  dominated  parabolic  type 
equations,  with  solutions  which  are  either  discontinuous,  or  with  discontinuous  deriva¬ 
tives,  or  containing  sharp  gradient  regions  which  are  difficult  to  be  completely  resolved 
on  todajPs  computer.  The  methods  we  investigate  fall  into  the  category  of  “shock 
capturing”  schemes,  which  means  that  these  methods  try  to  capture  shocks  or  other 
types  of  discontinuities  and/or  sharp  gradient  regions  with  a  relatively  coarse  grid, 
rather  than  completely  resolving  them.  These  methods  are  useful  when  it  is  either 
impossible  or  too  costly  to  completely  resolve  certain  solution  structure.  High  order 
accurate  finite  difference,  finite  element  and  spectral  methods  have  all  been  consid¬ 
ered.  Parallel  implementation  issues  of  these  high  order  methods  are  also  studied. 

Applications  of  the  algorithms  developed  and  analyzed  in  this  project  have  been 
carried  out  to  compressible  and  incompressible  flow  calculations,  and  semiconductor 
device  simulations.  We  have  been  using  the  Connection  Machine  CM-5  at  the  Army 
High  Performance  Computing  Research  Center  at  the  University  of  Minnesota  to  de¬ 
velop  highly  efficient  two  and  three  dimensional  flow  solvers  using  high  order  difference 
(ENO)  methods.  Collaborations  have  also  been  pursued  with  Dr.  Nisheeth  Patel  of 
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Ballistic  Research  Lab  and  Professor  Bernardo  Cockburn  of  University  of  Minnesota, 
about  joint  efforts  of  parallel  implementation  and  application  of  the  discontinuous 
finite  element  method  for  3D  Navier-Stokes  equations. 


3.  Summary  of  Research  Results 

Research  has  been  performed  in  all  areas  listed  in  the  original  proposal,  and 
progress  and  results  consistent  with  the  original  objectives  have  been  obtained.  There 
are  27  publications  (among  them  19  in  refereed  journals)  resulting  from  this  project, 
see  Section  4  for  a  list  of  them. 

The  main  class  of  finite  difference  methods  we  have  been  investigating  is  the  ENO 
(Essentially  Non-Oscillatory)  high  order  methods  using  fluxes  and  point  values,  first 
developed  by  Shu  and  Osher  in  1988,  using  the  adaptive  adaptive  interpolation  proce¬ 
dure  in  the  original  ENO  schemes  of  Harten,  Osher,  Engquist  and  Chakravarthy.  The 
advantage  of  the  version  of  ENO  schemes  based  on  point  values  and  nonlinearly  stable 
Runge-Kutta  time  discretizations,  designed  by  Shu  and  Osher,  is  that  the  resulting 
scheme  is  cost  efficient  for  two  and  three  dimensional  problems.  A  comparison  of  the 
efficiency  of  the  two  types  of  ENO  schemes  is  performed  in  [13]  (all  the  numbering  of 
references  are  according  to  that  of  Section  4).  In  this  same  paper  we  have  explored  the 
boundary  treatment  for  high  order  ENO  schemes  in  two  dimensional  problems  with 
a  non-rectangular  domain.  Both  wall  and  inflow/outflow  boundaries  are  considered. 
It  is  found  out  that  careful  treatment  of  different  boundary  conditions  is  essential  for 
the  realization  of  high  order  accuracy  and  stability  for  such  problems.  We  have  also 
found  out  that  the  grid  orientation  effect,  which  is  present  because  of  the  dimension 
by  dimension  feature  of  the  algorithm,  is  actually  diminished  with  the  increase  of  or¬ 
der  of  accuracy  and  is,  for  the  test  problems  we  have  computed,  negligible  for  fourth 
order  schemes. 

In  [11],  E  and  Shu  adapted  ENO  methodology  to  solve  incompressible  flows  and 
investigated  several  resolution  issues.  This  is  a  generalization  of  the  second  order 
Godunov  type  methods  used  by  Bell,  Colella  and  Glaz.  The  main  difference  from  the 
compressible  flow  is  the  presence  of  the  incompressibility  condition,  which  is  a  global 
condition  and,  if  not  treated  carefully,  may  affect  local  approximation  performance. 
Although  the  work  in  [11]  is  still  preliminary  (only  periodic  boundary  conditions 
were  considered,  for  which  the  incompressibility  condition  can  be  easily  enforced  with 
Fourier  transform  using  the  projection  method),  it  does  provide  the  key  evidence  that 
the  high  order  ENO  methodology  can  be  applied  to  incompressible  flow  and  has  good 
resolution  power  for  certain  important  physical  quantities. 

As  to  the  practical  issue  of  parallel  implementation  of  ENO  algorithm,  we  have 
successfully  developed  two  and  three  dimensional  ENO  code  for  the  Connection  Ma¬ 
chine  CM-5  at  the  Armj^  High  Performance  Computing  Research  Center  located  at 
University  of  Minnesota.  We  have  used  the  fact  that  communication  is  costly  and 
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other  tlia.a  next  neighbor  communication  should  be  minimized,  to  modify  the  im¬ 
plementation  of  the  algorithm  (without  changing  the  mathematical  content)  so  that 
communication  is  reduced  and  other  than  next  neighbor  communication  is  basically 
eliminated,  at  the  cost  of  a  slightly  increased  operation  count  and,  in  the  two  dimen¬ 
sional  case,  increased  storage  space  [3].  The  resulting  CM-5  code  using  128  nodes, 
for  a  two  dimensional  grid  size  of  400^,  is  almost  a  magnitude  faster  than  our  CRAY 
YMP  code.  Application  of  ENO  schemes  to  the  problem  of  two  dimensional  shock  in¬ 
teraction  with  a  hydrogen  bubble,  which  is  a  prototype  of  hypersonic  mixing  problem, 
is  also  carried  out  in  [3]. 

Adaptation  and  application  of  ENO  schemes  to  semiconductor  device  simulation 
has  been  actively  pursued.  There  are  various  models  in  device  simulation.  Most 
of  them  are  incomplete  parabolic  type  equations,  with  source  terms  coupled  to  a 
potential  equation,  and  with  a  significant  hyperbolic  part.  This  significant  hyperbolic 
part,  and  the  fact  that  solutions  can  typically  vary  by  several  magnitudes  over  a  short 
region  covered  by  a  only  a  few  grid  points,  warrants  the  application  of  high  order  shock 
capturing  methods.  Jointly  with  Professor  Joseph  Jerome  of  Northwestern  University 
who  is  an  applied  mathematician.  Professor  Umberto  Ravaioli  of  University  of  Illinois 
who  is  a  device  engineer,  and  their  collaborators,  we  have  investigated  two  dimensional 
energy  models  (both  the  hydrodynamic  model  and  the  energy  transport  model)  of 
device  simulation  from  both  a  modeling  and  a  computational  point  of  view.  My 
main  concern  in  this  project  is  to  understand  the  mathematical  structure  of  the 
various  models,  the  corresponding  solution  behavior  (especially  those  related  to  high 
gradient  regions  similar  to  compressible  gas  flow),  then  trying  to  resolve  issues  related 
to  the  adaptation  and  application  of  high  order  numerical  schemes  (both  ENO  and 
finite  element  techniques).  As  a  result  of  several  years’  research,  the  algorithm  based 
on  ENO  methodology  has  been  developed  into  a  robust  and  reliable  tool  for  device 
simulation  using  various  models.  In  [1]  a  new  energy  transport  model  is  discussed  and 
tested  to  show  its  improvement  upon  spurious  velocity  overshoot  at  the  right  junction 
for  a  standard  n-\-n-n+  channel;  in  [4],  [7]  and  [8]  the  hydrodynamic  model  and  the 
energy  transport  model  are  analyzed,  and  algorithm  based  on  ENO  methodology  is 
designed  and  tested  on  two  dimensional  problems.  This  algorithm  is  used  to  compare 
the  two  different  models  for  their  accuracy  in  describing  the  physical  phenomena  in  the 
device;  [18]  concerns  the  algorithm  development  using  discontinuous  Galerkin  method 
and  mixed  finite  element  method  for  the  hydrodynamic  model  in  two  dimensions,  and 
a  numerical  test  using  a  2D  MESFET  device;  [6]  studies  shock  formation  and  shock 
width  in  small  devices;  [19]  explores  the  response  of  the  hydrodynamic  model  to  heat 
conduction,  mobility  and  relaxation  expressions,  [15]  and  [24]  study  the  crucial  issue 
of  whether  transport  effect  is  important  or  not  in  one  and  two  dimensional  devices, 
the  result  is  positive  and  indicates  the  importance  of  using  the  full  hydrodynamic 
models  and  hyperbolic  based  algorithms. 

.Also  among  finite  difference  techniques  for  flow  calculations  is  the  class  of  compact 
schemes.  Compact  schemes  usually  can  be  designed  to  have  high  (spectral  like) 
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resolution  for  various  waves  and  are  particularly  suitable  for  long  time  calculation 
such  as  direct  numerical  simulation  of  turbulence.  In  [9],  Cockburn  and  Shu  provided 
a  framework  to  apply  compact  schemes  for  shocked  or  high  gradient  problems.  TVB 
compact  schemes  in  ID  and  L^o  stable  schemes  in  two  and  higher  dimensions  are 
discussed  and  tested.  The  basic  idea  of  compact  schemes  (using  rational  functions 
rather  than  polynomials  to  form  the  approximation)  can  also  be  combined  with  ENO 
ideas  to  obtain  another  type  of  ENO  schemes,  which  is  currently  under  investigation. 

The  combination  of  applied  analysis  and  direct  numerical  simulation  to  study 
physics  of  fluids  phenomena,  is  carried  out  jointly  with  Professor  Weinan  E  of  Courant 
Institute.  We  have  studied  the  effective  equations  and  the  inverse  cascade  of  energy  for 
the  two  dimensional  Kolmogorov  flow  in  [5],  and  the  small  scale  structure  in  Boussi- 
nesq  convection  in  [12].  Typically,  analysis  (formal  as  well  as  rigorous  asymptotic 
analysis)  is  carried  out  as  far  as  possible,  while  carefully  designed  numerical  simula¬ 
tion  (often  with  the  aid  of  partial  result  from  analysis)  is  carried  out  to  continue  the 
attack. 

The  important  issue  of  maintaining  positivity  of  pressure  and  density  in  high  or¬ 
der  finite  volume  type  compressible  flow  calculations  in  multidimensional  arbitrary 
triangulations  is  analyzed  in  [20].  Adaptation  of  ENO  schemes  to  solve  the  viscoelas¬ 
ticity  equation  with  fading  memory,  and  related  theoretical  study  of  shock  formation 
and  long  time  behavior,  are  performed  in  [26].  Application  of  ENO  ideas  to  computer 
vision  problems  is  performed  in  [27]. 

We  have  been  investigating  the  possibility  of  applying  spectral  method  to  shocked 
problems.  Two  approaches  are  considered;  The  first  is  of  non-oscillatory  type:  either 
to  add  into  the  spectral  basis  some  discontinuous  functions  to  make  the  approxima¬ 
tion  non-oscillatory  (Cai,  Gottlieb  and  Shu),  or  to  use  some  conservative,  high  order 
hybriding  between  spectral  and  ENO  fluxes  (Cai  and  Shu).  This  approach  has  the 
advantage  that  essentially  no  oscillations  will  appear  in  the  process  of  numerical  simu¬ 
lation,  hence  one  can  obtain  results  similar  to  those  obtained  by  high  resolution  finite 
difference  schemes.  Another  approach  is  to  use  some  mild  high  frequency  filtering  just 
to  stabilize  the  solution,  which  is  still  oscillatory,  then  to  post-process  the  solution 
to  obtain  accurate,  non-oscillatory  output.  This  approach  is  based  upon  the  phi- 
losoph}'  that,  although  the  numerical  solution  itself  may  be  oscillatory,  its  moments 
(i.e.,  its  integration  against  smooth  functions)  are  spectrally  accurate.  This  can  be 
proven  for  linear  problems  and  is  also  believed  to  be  true  in  some  sense  for  nonlinear 
problems.  In  order  to  perform  this  post-processing,  an  adequate  reconstruction  pro¬ 
cedure,  which  can  recover  spectrally  accurate  point  values  everywhere,  including  at 
the  discontinuity  itself,  from  the  oscillatory  spectral  partial  sum,  must  be  developed. 
Recently,  we  obtained  a  breakthrough  result  in  this  direction  [2],  [10],  [21],  [22],  [23], 
which  established  the  theoretical  foundation  of  this  approach  and  opened  the  road  for 
practical  computations.  The  recovering  of  exponentially  accurate  points  everywhere 
including  at  the  discontinuities  themselves,  is  rigorously  pro\en  and  successfully  im¬ 
plemented  even  for  small  to  moderate  N  (the  number  of  modes  or  grid  points),  for 
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botli  Fourier  and  polynomial  methods,  and  for  both  Galerkin  and  collocation.  The 
crucial  issue  of  assessing  information  content  of  using  a  spectral  method  to  solve 
a  nonlinear  conservation  law  is  studied  in  [25],  using  the  reconstruction  procedure 
described  before. 

Another  quite  successful  technique  for  hyperbolic  conservation  laws  is  the  discon¬ 
tinuous  Galerkin  finite  element  method.  In  this  method  the  partial  differential  equa¬ 
tion  is  multiplied  by  a  test  function,  integrated  over  a  cell,  and  formally  integrated  by 
parts  to  obtain  a  weak  formulation.  A  solution  is  sought  among  discontinuous  (across 
cell  interface)  piecewise  polynomials  of  r-th  degree  for  a  (r  +  l)-th  order  method.  Be¬ 
cause  of  the  discontinuity  at  cell  interface,  this  method  can  accommodate  successful 
finite  difference  methodology  (approximate  Riemann  solvers  and  limiters)  into  a  finite 
element  framework  (Cockburn  and  Shu).  Theoretical  results  similar  to  finite  differ¬ 
ence  methods,  such  as  total  variation  stability  for  ID  and  maximum  norm  stability 
for  2D  and  3D,  can  be  proved  for  this  class  of  discontinuous  Galerkin  methods  of 
arbitrary  order  of  accuracy  and  for  (almost)  arbitrary  triangulations.  An  essential 
difference  between  this  class  of  finite  element  method  and  the  finite  volume  method 
(which  can  also  be  defined  on  an  arbitrary  tri angulation)  is  that  the  latter  has  only 
one  independent  degree  of  freedom  (the  cell  average)  over  each  cell,  while  the  former 
has  many  (for  example,  it  has  three  degrees  of  freedom  for  the  piecewise  linear  case  in 
2D).  This  fact  renders  the  scheme  more  local  (no  wide  stencil  reconstruction  is  needed 
to  compute  the  flux  at  cell  interface),  hence  more  suitable  for  parallel  computing,  and 
provides  a  different  setting  for  theoretical  justification  of  stability  and  convergence  of 
the  algorithm.  For  example,  the  recent  result  of  Jiang  and  Shu  [14],  [16]  proves  a  cell 
entropy  inequality  for  the  square  entropy,  which  implies  L2  stability  in  general  and 
convergence  in  some  special  cases,  for  arbitrary  triangulations  in  any  space  dimension 
and  arbitrary  order  of  accuracy,  without  even  using  nonlinear  limiters.  This  result 
is  much  stronger  than  the  corresponding  results  for  finite  volume  or  finite  difference 
schemes,  which  typically  require  one  space  dimension,  a  convex  physical  flux  and  no 
more  than  second  order  accuracy  to  prove  the  same  entropy  inequality.  In  practice, 
finite  element  methods  can  handle  complicated  geometry  and  boundary  conditions 
more  easily. 
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